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Anderson localization of Bogoliubov excitations is studied for disordered lattice Bose gases in 
planar quasi-one-dimensional geometries. The inverse localization length is computed as function 
of energy by a numerical transfer-matrix scheme, for strips of different widths. These results are 
described accurately by analytical formulas based on a weak-disorder expansion of backscattering 
mean free paths. 


I. SETTING, OBJECTIVES, AND SCOPE 

Single-particle Anderson localization in quasi-ID ge¬ 
ometries with correlated disorder is already a challeng¬ 
ing problem [1- ]. With interactions, the situation be¬ 
comes even more interesting. Here, we study the localiza¬ 
tion properties of the Bogoliubov quasi-particles (BQP) 
of weakly disordered Bose gases at zero temperature on 
quasi-ID lattices. While interactions tend to screen the 
disorder and thus stabilize extended (quasi-)condensates 

] for weak disorder [6, 7], BQPs are expected to be lo¬ 
calized irrespective of their energy or disorder strength 
in low dimensions [8-10], thus emulating noninteract¬ 
ing particles in the orthogonal Wigner-Dyson universal¬ 
ity class [ ]. Yet, BQPs differ qualitatively from non¬ 
interacting particles because of their collective, phonon¬ 
like behavior at low energy. Moreover, they experience a 
randomness mediated by the inhomogeneous condensate 
background, which responds nonlinearly and nonlocally 
to the bare disorder [12-1- ]. This interplay has been ex¬ 
tensively examined in ID [15-18], where direct backscat¬ 
tering provides the only pathway to localization. 

In this Letter, we discuss the localization length of 
BQPs on quasi-ID planar lattices of transverse width 
N x 1 (shown in Fig. la). While the techniques devel¬ 
oped below also apply to N x x N y bars, strips provide the 
simplest realization of a multi-channel geometry, where 
phase-coherent scattering between modes is the rule. The 
bosons are described by the Bose-Hubbard (BH) model 
with on-site interaction U > 0, hopping t (with periodic 
boundary conditions across the strip) and on-site disor¬ 
der V x . The random potential V x is drawn from a Gaus¬ 
sian distribution and has no spatial correlation on the 
lattice scale, 1414' = 5 XX 'V 2 , where the bar denotes dis¬ 
order averaging and V 2 is the on-site variance. For each 
realization of disorder, the mean-held (MF) ground-state 
density n x is determined by the condensate amplitude 
<P X = \fnx that solves the discrete Gross-Pitaevskii (GP) 
equation 

-t X $x' + (2t + V x )$ x + U$ 3 x (1) 

( XX') 

where the sum runs over the nearest neighbors of x. The 
chemical potential fj, controls the average occupation n = 


nff. In the limit of large occupations [19], an expansion 
of the BH Hamiltonian around the MF solution yields 
the Bogoliubov-de Gennes (BdG) equations 
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where u x 

and v x are 

the BQP particle and hole 

com- 


ponents at energy E. These Bogoliubov excitations de¬ 
termine ground-state properties as well as the (thermo-) 
dynamical response of the system. 

In the remainder of this paper, we compute the de¬ 
pendence of the BQP localization length on the excita¬ 
tion energy. We first present a numerical transfer-matrix 
method that provides us with precise estimates of the 
localization length over a wide range of parameters (sec¬ 
tion II). Then we develop an analytical framework that 
combines random-matrix theory [1] with a microscopic 
transport theory (section III). In particular, we find that 
the predictions using a lattice Boltzmann mean-free path 
are in excellent agreement with the numerics. 

II. NUMERICAL TRANSFER-MATRIX 
CALCULATIONS 

Our numerical calculations of the localization length 
follow a standard transfer-matrix scheme [1, 20-22]: The 
BdG equations (2) are cast into the recursive form 
lOi+i = where the 4A-component vector Wi en¬ 

codes the values of u x and v x on two neighboring lattice 
slices at z = i and z = i — 1 (throughout the article, 
distances are expressed in units of lattice spacing). The 
transfer matrix M, is parametrized by E as well as the 
values of 14 and n x on slice i. The propagation of the 
BdG excitations along z, from an arbitrary initial con¬ 
dition Wq, then simply amounts to matrix-vector multi¬ 
plication. Via this procedure, the estimator ln||w£,||/.L 
yields, with probability one, the largest Lyapunov ex¬ 
ponent Ai in the limit L —> oo. The r largest Lya¬ 
punov exponents Ai > • • • > X r are retrieved similarly, 
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Figure 1. System configuration for the 4x1 geometry, (a) 
Sketch of the discrete strip geometry, with hopping ampli¬ 
tudes indicated on the left and on-site disorder on the right, 
(b) Clean lattice Bogoliubov dispersion showing contributions 
from the transverse modes kj, (j = 0,... ,3). (c) Clean ex¬ 
citation DOS [cf. Eq. (3)] with ID divergencies at channel 
openings. The color code indicates the relative contribution 
of the individual transverse modes. All distances are given 
in units of lattice spacing. Panels (b) and (c) are drawn for 
Un/t = 2.0. 


by propagating a frame of r linearly independent initial 
vectors w\^ (j = 1 ,...,r), enforcing their orthogonal¬ 
ity during propagation, and monitoring their exponential 
growth (to maintain numerical accuracy, we performed a 
Gram-Schmidt orthonormalization every 8 sites). As the 
BQP Lyapunov exponents come in pairs of opposite sign, 
2TV vectors are required to calculate the smallest positive 
Lyapunov exponent A = X- 2 N, identified as the inverse lo¬ 
calization length 1/Zioc [1, 22]. 

One crucial advantage of the transfer-matrix approach 
is the self-averaging of the In ||u?j^||/L estimators, whose 
relative fluctuations decrease as 1/ y/L. The propagation 
of one set of initial vectors over a large distance L thus 
suffices to compute the Lyapunov spectrum with the de¬ 
sired precision. Such a scheme is easy to implement with 
an uncorrelated random potential, as the purely local 
random values of V x can be generated on the fly. In the 
case of BQPs, however, the GP equation (1) implies a 
global minimization problem that is already intractable 
for the size required to achieve only 10% accuracy on the 
Lyapunov exponents. In our calculations, we chained 
quasi-ID segments of length L 0 = 2 14 and, on each seg¬ 
ment, solved Eq. (1) at fixed density n with a conjugate- 
gradient technique [23, 24]. While each segment inter¬ 
face introduces a slight mismatch in //, and also in n x 
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Figure 2. Inverse localization length A (error bars indicate 
the error of the mean after averaging over 500 configurations 
in (a) and (b) and 225 in (c), and are hardly appreciable on 
the scale of the plot) and its estimates from transport theory 
as a function of energy. The dashed line is [(Ne + l)^bs ] _1 
from Eq. (5), and the solid red line is [2 (Ne + l)hr ] _1 from 
Eq. (6). The staircase-like line indicates the number of open 
channels Ne (scale on the right border). Parameters: hopping 
Un = 2t ~ /i, i.e., the transition from sound waves to particle 
excitations is around E/t = 2. Disorder strength V = 0.1 Un. 


over typically yjt/(Un ) sites, the impact of these arti¬ 
facts on A vanishes with increasing Lq. We checked that 
all the results presented here are converged with respect 
to the choice of Lq. To ensure a sufficient scrambling 
of initial data and speed up estimations, we chained up 
to 64 segments and averaged the In ||m^|| data over 200 
independent propagations. 

Figure 2 shows the inverse localization length A as func¬ 
tion of energy for various strip widths N. For all N we 
observe two important features. First, A diverges at the 
inner band edges, where the channels of the clean sys¬ 
tem open or close (cf. Fig. lc) and thus contribute a 
divergent ID density of states (DOS), which results in 
strong scattering and thus localization on short length 
scales. Second, for an equal number of open channels 
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and once inner divergences are disregarded, A displays 
a similar trend as in ID [10], increasing with E in the 
phonon regime E < Un, and decreasing in the free- 
particle regime E >Un. 


III. ANALYTICAL TRANSPORT THEORY 


For quasi-ID strips, random-matrix theory [1] predicts 
that the localization length Zi oc = 2 (Ne + l)Ztr is propor¬ 
tional to the transport mean free path l tr . Here, Ne is 
the number of open channels at energy E. In this section, 
we validate this prediction using two different analytical 
expressions for the mean free path. 

The clean system is discrete translation invariant, and 
therefore the free-particle dispersion relation is diagonal 
in momentum, = 2t£A_ a , z [l — cos (ft*)]. The dis¬ 
persion of Bogoliubov excitations then reads as usual 

= [ 4(4 + 2 Un)] 1 / 2 . On the strip, the momentum 
k = (kj } k z ) is transverse quantized as kj = 2irj/N un¬ 
der periodic boundary conditions with j = 0,..., N — 1, 
thus defining the N possible channels. In terms of the 
longitudinal group velocity Vk = dek/dk z (all velocities 
here and below refer to the ^-component), the DOS per 
unit length and transverse size is expressed as 


p{E) = 



dk z 
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( 3 ) 

Se denotes the energy shell in momentum space, namely 
the set of points k such that £& = E. Its cardinality 
\Se\ = 2 Ne counts both forward and backward propa¬ 
gating channels with k z > 0 and k z < 0 , respectively. 
The clean DOS is plotted in Fig. lc. The opening and 
closing of each channel produces a characteristic van 
Hove divergence. Near the divergence, essentially only 
one channel matters. In between these singularities sev¬ 
eral channels mix, as indicated by the color coding in the 
plot. 

The disorder potential produces elastic scattering out 
of a given mode k with a rate 7 k = ^Ylp^( e k — 


£p)\Wk P \ 2 . Here W kp = V k - p w k p is a (first-order) 
matrix element of the disordered Bogoliubov Hamilto¬ 
nian, namely the bare disorder Fourier component Vk- P 
dressed by an envelope 


Wkp — 


(4 + Un Vk- P ~ 2Une° k 
e k (2Un + 4 _ p ) 


(4 = 4 ) ( 4 ) 


that accounts for the underlying nonlinearities [13, 21]. 
[Eq. (4) is the on-shcll value | e °= e ° given by Eqs. 
(13) - (15) in [25] or, equivalently, the lattice version 
of Eq. (14) in [26].] The associated elastic mean free 
path is then the group velocity divided by the scattering 
rate. Since localization along the 2 direction can only 
be caused by scattering events that flip the momentum 


^-component, we consider the backscattering mean free 
path 


WEr ‘ = s - < 5 > 

The primed sum runs only over those momenta p £ S Ck 
with z-conrponents of opposite sign, sgn (p z k z ) = — 1 . 
Moreover, we take the shortest inverse mean free path 
at a given energy E since the numerics produce the 
smallest localization exponent. The resulting predic¬ 
tion [(Ne + l)Zbs(-E ')] -1 is plotted as the dashed line in 
Fig. 2. There is already some overall qualitative agree¬ 
ment, which becomes even quantitative at low and high 
energies, where only a single channel contributes. And 
indeed, in this case (Ne = 1) one has 2nNp(E) = 2\vk\~ x 
at £fc = E, and thus the estimate for the Lyapunov ex¬ 
ponent ( 2 lbs )” 1 = 1 i / 2 ' u, fe(-fc)/( 2 A r u^) agrees with the es¬ 
tablished theoretical result [10, 13, 14]. 1 

The number Ne of open channels is also plotted in each 
figure. With more open channels, the estimate involving 
lb s becomes less reliable, which becomes especially obvi¬ 
ous near the center of plots b) and c) for the 8 x 1 and 
32 x 1 geometries. The agreement deteriorates even fur¬ 
ther if only diagonal backscattering into the same channel 
(pj = kj and p z = — k z ) is taken into account. 

We therefore try to improve on our estimate and con¬ 
sider the lattice version of the full Boltzmann transport 
mean free path, which can be derived within linear- 
response quantum transport theory [27, 28]. With the 
notations introduced above, this gives 

Vk)w 2 kp 
Vk\ 

( 6 ) 

Here, (v 2 ) E = K|/[27rY p(E)\ denotes an 

energy-shell average. The usual factor (1 — cos(?) in 
scattering angle 6 under the momentum integral (see, 
e.g., [13, Fig. 6 ] and [29, Eq. ( 6 )]) is now represented by 
(v p — Vk)- The factors of and \v p \ in the denom¬ 
inator stem from the A^-integrals over on-shell spectral 
functions, analogous to the second equality in (3). 

For a strict ID problem, i.e., for a single open chan¬ 
nel (Ne = 1 ), this expression reduces to 2 l tI = lbs, as 
it should [1, Eq. (152)]. But also when several chan¬ 
nels are open, the agreement with the numerical data is 
clearly excellent, as shown in Fig. 2. Remarkably, Eq. ( 6 ) 
performs especially well in between the ID resonances, 
where the coherent coupling between channels has to be 
described rather accurately. 


ltr(E)~ = 


V 2 


2ttN 2 p(E)(v 2 ) 3 J 2 


E 


U P (^p 


E k,p£SE 


1 In the ID case with zbfc along z and quadratic dispersion 
e 2 k = 4e fc’ one has w k{-k) = e fc/ e fc and thus [2Z bs (e fe )] _1 = 
V 2 Sl/[2N(vl) 2 } in terms of the free-particle velocity v = 
de^/dk and Sk = e®/(£/n + €°), which agrees with [10, Eq. (11)]. 
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Of course, the result (6) is perturbative in nature, es¬ 
sentially valid for weak enough disorder (V <C Un ), and 
thus only holds if the (quasi-)condensate is extended. 
For much stronger disorder, the condensate fragments 
and the Bose fluid enters the non-superfluid, Bose-glass 
phase [24, 3i ]. 


IV. SUMMARY AND OUTLOOK 

In conclusion, we have calculated the localization 
length for Bogoliubov excitations of Bose-Einstein con¬ 
densates on disordered lattices with a planar quasi-ID 
geometry. Numerical data from a transfer-matrix com¬ 
putation are very well reproduced by an analytical for¬ 
mula adapted from continuum transport theory. Our ap¬ 
proach could extend to larger systems N x 1 (and N x N) 
with N —> oo using finite-size scaling and thus yield the 
full 2D (and 3D) localization lengths, with the aim to 
track down possible mobility edges for BQPs in dimen¬ 
sions higher than 1. 
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